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Motivated by some inconsistencies in the way quantum fluctuations are included beyond the 
classical treatment of hard-core bosons on a lattice in the recent literature, we revisit the large-S 
semi-classical approach to hard-core bosons on the square lattice at T = 0. First of all, we show 
that, if one stays at the purely harmonic level, the only correct way to get the 1/S correction 
to the density is to extract it from the derivative of the ground state energy with respect to the 
chemical potential, and that to extract it from a calculation of the ground state expectation value 
of the particle number operator, it is necessary to include 1 / \/~S corrections to the harmonic ground 
state. Building on this alternative approach to get 1 / S corrections, we provide the first semiclassical 
derivation of the momentum distribution, and we revisit the calculation of the condensate density. 
The results of these as well as other physically relevant quantities such as the superfluid density 
are systematically compared to quantum Monte Carlo simulations. This comparison shows that 
the logarithmic corrections in the dilute Bose gas limit are only captured by the semi-classical 
approach if the 1/5* corrections are properly calculated, and that the semi-classical approach is able 
to reproduce the 1/fc divergence of the momentum distribution at fe = 0. Finally, the effect of 1/5'^ 
corrections is briefly discussed. 



PACS numbers: 

I. INTRODUCTION 

Models of interacting bosons on a lattice are ubiqui- 
tous. They have been introduced to describe the low 
energy physics of systems as different as thin supercon- 
ducting films [1], Josephson junction arrays [2], '*He on 
substrates [3-5], cold atoms in optical lattices [6], bipo- 
larons [7], or quantum magnets in a field [8, 9]. In simple 
(unfrustrated) geometries quantum Monte Carlo (QMC) 
simulations do not suffer from any minus sign problem, 
and the resulting picture is often quite clear [10]. How- 
ever, in many recent applications, the relevant effective 
model contains terms that lead to a severe minus sign 
problem. This is for example true for the bosonic de- 
scription of frustrated quantum magnets in a magnetic 
field [11], where QMC approaches are not appropriate. 
To investigate such models, it is important to develop 
alternative approaches. 

A very important subclass is that of models of hard- 
core bosons on a lattice in which the on-site repulsion 
is assumed to be infinite so that it is impossible to have 
more than one boson at a given site. Such models appear 
for instance very naturally in the description of dimer- 
based spin- 1/2 quantum magnets in a field [9]. In this 
paper, we will concentrate on a model of hard-core bosons 
on a two-dimensional square lattice described by the sim- 
ple Hamiltonian: 

"H = -t ^(a|aj + ajfl]) - /^^n^, (1.1) 

where t is the hopping amplitude between neighboring 
sites, n is the chemical potential and a| (a^) denotes the 
operator creating (destroying) a hard-core boson at site 



i. One distinct property of hard-core boson models is 
that they can be mapped exactly onto spin-1/2 models 
using the Matsuda-Matsubara transformation [3]: rii = 
Sf + 1/2, al = S+ and a, = S'r. The equivalent spin-1/2 
model is a ferromagnetic XY model with magnetic field 
fi pointing in the z direction: 

n^-tY,2 [S^S] + sfs]) is^ + 1/2) (1.2) 

(id) 

On the basis of this mapping, a semi-classical approxima- 
tion can be developed starting from the large S limit of 
this spin Hamiltonian. This approach has been developed 
in a series of papers [12-14]. The paper by Bernardet 
and coworkers [14] includes a careful comparison with 
QMC simulations and shows that, already at the order 
of linear-spin wave theory, the semi-classical approach is 
quantitatively accurate. 

Building on this success, this semi-classical approach 
has recently been used quite systematically in the in- 
vestigation of frustrated models [13, 15-20] for which 
it is often the only available analytical approximation. 
These studies have revealed a number of subtleties how- 
ever in the implementation of the semi-classical approx- 
imation. A recurrent problem concerns the calculation 
of the bosonic density as a function of the chemical po- 
tential [21, 22], or equivalently of the magnetization as 
a function of the field. Bernardet et al. have calculated 
the density as the opposite of the derivative of the energy 
with respect to the chemical potential, which is equiva- 
lent to calculating the magnetization as the opposite of 
the derivative of the energy with respect to the field. But 
one could in principle equally well calculate the magne- 
tization as the expectation value of the operator in 
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the ground state. However, at the harmonic level, the 
two definitions do not lead to the same answer, and it 
is not clear which definition should be preferred. In ad- 
dition, in its current setting, the semi-classical approach 
only allows to calculate in a systematic way quantities 
that can be derived from the ground state energy, i.e. 
the density and the superfluid stiffness. For instance, no 
attempt has been made so far to calculate other ground 
state properties such as the momentum distribution func- 
tion. Finally, no attempt to check the convergence of the 
1/5 expansion by calculating higher order corrections has 
been made. 

In the present paper, we address all these issues. First 
of all, we show that, to get the same result using the two 
definitions of the density, one has to include l/VS cor- 
rections to the harmonic ground state. These corrections 
have to be included to get the correct answer to order 
l/S because they contribute at this order when calculat- 
ing the expectation value of S^. In the low density limit, 
we also show that these corrections are crucial to get the 
logarithmic corrections predicted long ago for interact- 
ing two-dimensional bosons. Building on this success, we 
use this corrected ground state to calculate the momen- 
tum distribution function, and we show that it leads to 
a divergence at zero momentum that agrees with QMC 
results. We also provide two complementary ways to cal- 
culate the condensate at the order 1/5: from the deriva- 
tive of the energy with respect to a transverse field, and 
by a calculation to the zero-momentum occupation factor 
using the perturbed ground-state. Finally, we calculate 
the 1/5^ correction to the ground state energy and shows 
that it improves over the 1/5 result, supporting the ba- 
sic assumption of the semiclassical approach that the 1/5 
expansion is well behaved even for 5 = 1/2. 

Let us emphasize that we agree with all the results of 
Rcf. [14] to order 1/5. In that respect, the main objec- 
tive of the present paper is to show how these results 
can be obtained from perturbing the harmonic ground 
state, with two new results: a clear answer regarding 
the appropriate way to calculate the expectation value of 
observables at the order 1/5, and the first semiclassical 
calculation of the momentum distribution. 

This paper is organized as follows. In Sec. II the 
model is treated in the context of linear spin wave theory. 
Sec. Ill is devoted to the semiclassical correction of the 
harmonic ground state, and to the computation of sev- 
eral observables in this perturbed ground-state. Sec. IV 
presents a comparison of the spin wave results obtained 
with the results of QMC simulations. Sec. V discusses the 
validity in the context of a 1/5 expansion of the sum rule 
which states that the total density is equal to the sum 
of the condensate density and of the average momentum 
distribution function. Sec. VI presents some results ob- 
tained beyond the linear spin wave approximation. A 
short conclusion is given in Sec. VII. Finally, some details 
about the calculation of the superfluid density and of the 
momentum distribution are given in Appendices A and 
B. 



II. LINEAR SPIN WAVE THEORY 

A. The model 

To perform a semi-classical expansion, it will prove use- 
ful to extend the model of Eq.(1.2) in two ways. First of 
all, we rescale the amplitudes in such a way that the vari- 
ous terms are of the same order in the large 5 limit while 
the Hamiltonian of Eq.(1.2) is recovered for 5 = 1/2. 
Secondly, and more importantly, we introduce a trans- 
verse field r > in the x direction. These modifications 
lead to the Hamiltonian: 

11 _ _ ^ ( qxqx , qV qy\ \^ qZ ^ 3. 

i i 

(2.1) 

The introduction of a transverse field F turned out to 
be an essential ingredient in two respects. On one hand, 
it allows one to calculate the condensate density as the 
opposite of the derivative of the ground state energy with 
respect to F, hence to get an expression that is correct to 
order 1/5. On the other hand, it breaks the continuous 
U{1) symmetry of the Hamiltonian of Eq. (1.2) and opens 
a gap in the spectrum of the model. Thanks to this 
gap, the correction to the harmonic ground state is not 
divergent, and the corrected ground state can be used to 
calculate the expectation value of various observables to 
order 1/5. The results for the original model are then 
obtained by taking the limit F — > of the expectation 
values. 



B. Classical solution 

In the classical limit, spin operators are replaced by 
three-dimensional vectors of norm 5. In the absence of a 
transverse field, the ground state consists of spins ordered 
ferromagnetically in the x — y plane with a longitudinal 
magnetization m that varies linearly with the magnetic 
field /Lt until saturation. When F > 0, the classical solu- 
tion lies in the x ~ z plane and can be parametrized as 
follows: 







/ sin 


sf 











I cos 



(2.2) 



With this parametrization, the classical energy per site 
is given by: 

= -2tsin2 6l-/icos6l-Fsin6l. (2.3) 

The angle 6 is fixed by minimizing the classical energy, 

- At sin e cos + fi sin 6* - F cos 6* = 0. (2.4) 

In the limit F — > 0, Eq. (2.4) has the simple solution 
cos^o — m/4^, and Sf is different from zero in the field 
range —4 < /i/t < 4. This defines the critical magnetic 
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field /ic = —4:t at which the system starts to acquire a 
transverse magnetization. For S = 1/2, at fi > /Zc, we 
have S'f > —1/2 which in terms of the original hardcore 
boson model corresponds to a non zero density of bosons. 
Thus fic is the value of chemical potential at which the 
system is no longer empty and an hardcore boson popula- 
tion starts to develop. For F 7^ 0, the angle 6* is a function 
of F. In the following, we will focus on small F case, and 
we will calculate the small F correction to several quan- 
tities. From the equation sin 0(— 4t cos 9 + fJ.) = F cos 9, it 
is easy to see that the first order correction to 6 is given 
by: 



r=o 



cos 6*0 
U sin^ 6*0 ' 



(2.5) 



The classical on site magnetization is — Scosd and 
the classical hardcore boson density is given by p"^^^^^- 
(cos 9 + l)/2. In the limit where the transverse field van- 
ishes the classical density is given by: 



^class. 



if^ + 1 

2 \4t 



(2.6) 



C. Holstein-Primakoff transformation 

In order to study the effect of quantum fluctuations 
around this classical solution, we start by performing a 
rotation of the spins at each site 



Sf = cos 9Sf - 

CV _ QV' 

Sf = - sin 9Sf 



sm9Sf 
+ cos 9Sf 



(2.7) 



such that the Hamiltonian of Eq. (2.1), expressed in 
the rotated frame (x', y' , z'), has a ferromagnetic ground 
state. The new spin operators can be expressed in terms 
of Holstein-Primakoff bosons [23]. To next to leading 
order, the expressions take the form: 



Sf 
Sf 

sf 



= S-b]b, 



'2S_ 
"2~ 
2S 



2i 



1 



4V2S 
1 



4:W2S 



nib. 



bjn. 



b]n, 



(2.8) 

The resulting Hamiltonian in terms of Holstein-Primakoff 
bosons can be expanded as 



(2.9) 



where "H^"-* is proportional to 5 3 . The first term of this 
series is H'^^^ = A^i?(o), N being the total number of sites. 
By construction, 'H'^-' = since we expand around a spin 
configuration which is a classical minimum of the energy. 



■H^^^ is quadratic in bosonic operators while "H^^-* and 
H^^^ contain only three or four boson terms respectively. 
Their expressions are given by: 



2S 



Y,icos'e~l){b,b^+blb^^ (2.10) 



^^blb, {4t sin^ 6* + cos 6* -I- F sin 61) , 



^(3) 



^(4) ^ _ 



2t 



SV2S 



J^niibj +b])sm9cos9, (2.11) 



(hi) 



t 

S^ 



^ ^(1 - cos^ 9) {[n, + nj]bibj + h.c.) 
5]-^(l + cos2 0) (b\[n,+n,]b^+h. 
l^^^sin^ 9 niUj. (2.12) 



t 
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^2 



The calculation of 1/ S corrections, to which most of the 
paper is devoted, is based on T^'^) ^nd 'H'^^\ The fourth 
order correction H^^' will only be used in Section VI when 
we calculate the l/S"^ correction to the energy. 



D. Diagonalization of the harmonic Hamiltonian 

In terms of the Fourier transformations of the Holstein- 
Primakoff operators defined by 

b, = ^y foke^^^" b\^^y ble-^^^"^ 



N 



(2.13) 

'H'-2) can be decoupled into a sum over different modes. 
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2S 
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^ ^(4tsin2 9 + /i cos 6* + F sin ^ 



(2.14) 



where the coefficients ^k and Sk are defined by: 
Ak = -^7k(cos^6' + 1) + 2isin2 6* + ^ cos6' + ^ sin( 



t 2 
Bk = ^DiSm 9, 



(2.15) 



with 7k = cosfca; -I- cos fey. With the help of Eq.(2.5), 
these coefficients can easily be expanded to linear order 
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in T: 



and 



r 

2 



cos^ 6>o 7k ^ 1 
sin 00 2 sin 6q 



(2.16) 



4 sin 6*0 



where 



[7k (l + cos2 6*0) - 4] /2 and B£ = 
[t7ksin^ ^o] /2 denote the coefficients and Bk in the 
absence of a transverse field [14]. The second term in 
Eq. (2.14) can also be expanded to first order in F, lead- 
ing to -(l/S) Ek 2i + r/(2sin0o)]. 

The quadratic Hamiltonian (2.14) can be diagonalized 
via a Bogoliubov transformation: 



bk Ukak - ^'kolk 



Wk^-k- 



(2.17) 



The coefficients which ensure that the operators «k('^k) 
satisfy bosonic commutation relations and that the 
Hamiltonian is diagonal are given by: 




(2.18) 



In terms of the Bogoliubov operators, and to first order 
in r, the Hamiltonian takes the diagonal form: 



-T 

k 

-y 



2t 



2 sin e. 



(2.19) 



The ground state of the harmonic Hamiltonian "H*-^-* is 
the vacuum of a particles. We will refer to it as the 
harmonic ground state in the rest of this paper. The 
first order l/S correction to the energy per site is given 
by: 



^(2) = J_y 

SN ^ 



^ ^ 2 sm ( 



(2.20) 



The only difference with the approach of Ref. 14 is that, 
as long as the transverse field is strictly positive, the Bo- 
goliubov transformation is well behaved even at k = 
since the excitation spectrum is gapped. Indeed, for small 
k and F, the excitation energy 51k = 2^y4^ — Bl/ S can 
be written as 



(2.21) 



with 



A= -\/nsin0o + O(r5), 



V = — isin^o 



r(l + 3cos2 6 



©(r^). (2.23) 



In the limit F — >■ 0, the spectrum becomes gapless and 
linear, as expected for phonon-like excitations in a super- 
fiuid. 



E. Calculation of the densities from the ground 
state energy 

A system of bosons is characterized by three densities: 
the total density, the condensate density, and the super- 
fiuid density. They can all be calculated as derivatives 
of the ground state energy. Using the Hellman-Feynman 
theorem which states that 



(dn{h)\^ 



(2.24) 



where h is some parameter of the Hamiltonian, one can 
calculate the longitudinal magnetization m as 



m{S) ^ 



a£;(2)(r = o) 

9/i 



and the transverse magnetization mj_ as 



m_L(5) 



-S 



a£;(2) 



dV 



(2.25) 



(2.26) 



r=o 



while the spin stiffness is given by the second derivative of 
the energy with respect to a twist (see Appendix A). The 
advantage of deriving these densities from the ground 
state energy is that, once we have an expression of the 
energy to a given order in l/S*, we obtain expressions of 
the densities which are correct at the same order. Let us 
discuss the result for the various densities. 



1. Total density 

Using the expression of the energy of Eq.(2.20) for 
F = 0, the derivative with respect to jjL leads to the lon- 
gitudinal magnetization 



The total density p is related to the longitudinal magne- 
tization hy p = m{S = 1/2) + 1/2, which leads to 



1 + COS0O , COS0O 1 Mk - B'^ 



4 N ■ 



K + B^ 



(2.22) perfect agreement with Ref. 14. 
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2.. Condensate density 

Taking now the derivative of the energy of Eq.(2.20) 
with respect to F, we obtain the following expression for 
the transverse magnetization: 



fluctuations (see Eq. (2.20)), and the second one is due 
to the increase of the total density of particle p as seen in 
Eq. (2.27). Interestingly, the same term appears in the 
condensate density, but with an opposite sign. 



2sin6'o N ^ 

k 



- 1 



cos (70 1 ■^-^ 

4sin6'o^^ 

k 



7k 1 



The condensate density which is the number of bosons 
occupying the k = mode per site po = (Q^k=o'^k=o) 1^ ~ 
(S^ Sj) /N^ , is simply related to the transverse 
magnetization by po — [m±{S = 1/2)]^, which leads to 
the expression 



Po 



sm Oo ^ 



1 

2N 



E 



- 1 



cos^ 9o 1 

4 N 



E 



7k A 



B'u 



B^- 



(2.28) 



This expression is different from that of Ref. 14. The two 
expressions are strictly equivalent only at p = —At (low 
density limit) and at p — (half filling). In the range 
—4<p/t<0 and < p/t < 4 the expression of Ref. 14 
differs form that of Eq. (2.28) by a term which is of order 
l/S*^ [33]. We believe that the above expression for the 
transverse magnetization is the correct one to order 1/5. 
This will be further supported by a direct calculation of 
the expectation value of S'^ in the next section. 



3. Superfluid density 

As explained in the appendix A, the superfluid density 
is given at 1/S order by the following expression 



Psf 




2i-v/(^£)^-(^k)') 



B^- 



(2.29) 



This expression is equivalent to order 1/5 to the expres- 
sion of Ref. 14. Looking at the expressions Eqs. (2.28) 
and (2.29) for the condensed and superfluid densities, one 
can make several interesting observations. First, quan- 
tum fluctuations deplete the condensate whereas they en- 
hance superfluidity. In Eq. (2.29), two contributions in 
the mechanism of superfluidity enhancement are present: 
the flrst term comes from the nearest neighbor kinetic en- 
ergy which increases (in absolute value) due to quantum 



III. LARGE S CORRECTIONS TO THE 
HARMONIC GROUND STATE 

Now that we have expressions for the longitudinal and 
transverse magnetizations valid to order 1/5, let us show 
how these expressions can be obtained as expectation val- 
ues of 5^ and 5^. Since the expressions to order 1/5 have 
been derived from the energy calculated at the harmonic 
level, one might expect that it is sufficient to calculate the 
expectation value of 5^ and 5^ in the harmonic ground 
state. As we shall see, this is not the case. The basic 
reason is quite simple: in terms of Holstein-Primakoff 
bosons, the operators 5^ and 5^ contains terms of or- 
der 0(5) and terms of order 0(1). To get an expression 
which is correct up to order 0(1), i.e. which includes 
all corrections up to 1/5, one should thus include in the 
ground state corrections up to order 1 / \/5, if any, since 
the expectation value of the 0{S) part of the operators in 
such a correction will give a contribution of order 0(1). 
As we shall now show, the term H^^^ in the expansion 
of the Hamiltonian indeed leads to a correction to the 
ground state of order 1 / y/S. 



A. Beyond the harmonic ground state 

The objective of this part is to compute the large 5 
corrections to the harmonic ground state. To do so, we 
treat H^^^ Eq. (2.11) as a perturbation to 'H^'^\ the smaU 
parameter being 1/5. The ground state of H^^-*, the vac- 
uum of a quasiparticles, being non degenerate, we use 
Rayleigh-Schrodinger non degenerate perturbation the- 
ory. To first order in perturbation, the perturbed ground 
state is given by: 



10) +E 



|0) 



El 



|0) 



Si 



{e\H 



(3)1 



0{S) 



0( 



sVs' 



(3.1) 



where |0) denotes the vacuum of a quasiparticles, |e) 
magnon excitations and E^q^ iE\e)) the energy of the vac- 
uum (excited states). The first correction to the ground 
state is of order 1/y/S since l/(i?|o) — E^^-^) is of order 
0(5) while (e|H(3)|o)|e) is of order 0(l/5i). The sec- 
ond fine defines It/f^), the ket that gives the l/VS cor- 
rection to the ground state. The above wave function is 
correct to order l/\/5 since all terms "H*^"^ with n > 4 
in the Holstein-Primakoff expansion (2.9) will contribute 
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corrections of higher order in l/S. Finally, we have to 
normalize the state, which leads to: 

where C = ((/)3|05) and (V'olV'o) = 1 + 0{1/S'^). To 
compute \(f>^), we first express 'H'-'^-' in Fourier space 

(3.3) 

with 

^q+k^k^q = (^k+q4+q - "k+qa_k-q) (3-4) 

X (wkak - '"kalk)("q"q ~ ■^q"-q) 



(S^) in terms of the spin operators in the rotated frame, 
and to use the expansion of these operators in terms of 
Holstein-Primakoff bosons. This leads to: 

(Sf) = cos0{Sf)~sine{Sf) 

= cos 6 (^S - {blb,)^ 

^sinei^ih + bl)- ^{nA + bU)) 

(3.8) 

At the classical level, the average magnetization is given 
by Scos9o when F = 0. The spin wave corrections to 
this result are of order 0(1). Thus, given the structure 
of the perturbed ground-state IV'o) « (1 - C/{2S)) |0) + 
5~2|02)j the terms entering the average magnetization 
to order 0(1) are: 



4-q^k&q = K-q4-q - ^^k-qa_k+q) (3-5) 

X (wk^k - '^kalk)("q"q ~ ^qO^-q)- 

"H*^^^ being a three-body operator, its effect on the vac- 
uum is to create one-magnon or three-magnon excita- 
tions. In Sec. IIIB, we will show that, for the computa- 
tion of first order corrections to the average values of the 
observables of interest in this paper, only single magnon 
excitations are relevant. We thus write ji^s) as a sum of 
one-magnon and three-magnon contributions: 



\4>^) = |05)lm + |0')3m 



(3.6) 



and we concentrate on the expression of the one-magnon 
contribution \4>^)im- Due to momentum conservation 
in Eq. (3.3), the only single particle excitations allowed 
have zero momenta. Hence, E^^y^ — _E|g) = — 

-'iyjAl - Bl/S and \(j)^)im is given by: 



I Im 



V2 



sin 9 cos 9 



E 

k 



2vt 



{uq - Vq) 

- 7k (wk 



WkUk)] |lq= 



(3.7) 

where |lq=o) denotes an excited state of one magnon 
with momenta q = 0. Note that it is only possible to 
write down such an expression because we have included 
a transverse field in the Hamiltonian, so that the Bo- 
goliubov transformation is not singular at k = 0. This 
expression actually diverges in the limit F — ^ because 
y^Al - Bl = 0(Fi/2) while {uq - vq) = C'(Fi/4). As we 
shall see, the limit F — > must be taken after calculating 
the expectation value of the operators. 



B. Expectation values of observables 

1. Total density 

Let US first use the perturbed ground state to calculate 
the expectation value of S^. The first step is to express 



{Sf)\r=Q = 5cos6io - lini |cos6i(0|6|5,|0) 

+ sin0^((O|6l-|-&,|05)i„+h.c.)}. 

(3.9) 

In the above expression, we have only included \(l)^)im 
in the matrix element of bl + bi since the operator 
bl + bi can at most create or destroy one Bogoliubov 
excitation. The three magnon component \(t)^)3rn would 
only contribute to the matrix element {l/^/S){0\{nibi + 
bln,)/{4,V2S)\4>^), but this term is of order 0(1/5*) and 
can be neglected since we are interested in the 0(1) cor- 
rection to the expectation value of 5^ . 

The matrix elements (0|6|6J0) and {0\bl + b^\(j)^)im 
are readily computed in Fourier space: 



(3.10) 



and 



N 



{uo - VQ)aa\4)^)im 



t 

71 



smf cost/— ^ 



1 -f 7k (w^ - WkWk) 

N 



with 



Aa + Bo = 2t sin^ 9o + T 



cos^ 9o 



Ao + Bo 



1 



(3.11) 



0(F2). 



sin do 2 sin 6o 

(3.12) 

Note that the expression of the second matrix element 
has a finite F — > limit because of the extra uq — vq fac- 
tor. Injecting back Eqs. (3.10) and (3.11) into Eqs. (3.9), 
one recovers exactly the expression of the longitudinal 
magnetization obtained before from the derivative of the 
energy calculated at the harmonic level. 
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2. Condensate density 

The same procedure can be repeated for the operator 
(5f) = sin 6'(S'f') + cos9{Sf), and the T ^ hmit of its 
expectation value is given by: 

(5f)|r=o = 5sin6lo- hni {sin6i(0|6l6jO) 

-cos0^((O|6l+6J05)i„ + h.c.)}. 

(3.13) 

Injecting back Eqs. (3.10) and (3.11) into this expression 
leads to exactly the same expression for the condensate 
as the one obtained from the derivative of the energy with 
respect to F. 



3. Momentum distribution 

The main advantage of this approach is that it gives 
access to observables that cannot be calculated as deriva- 
tives of the energy. Among them, a physically very im- 
portant one is the momentum distribution defined by: 



(<«k) 



hm — 
r^o N ^ 



{sts-) 



(3.14) 



The details of the calculation of {S^Sj ) in the perturbed 
ground state (3.2) are given in Appendix B. For k ^ 0, 
the momentum distribution is given by: 



(4«k> 



1 , 

-(1 + cos( 
4 



cos^ 0o)vl 



ui^vi^ sm 



.(3.15) 



The classical expression for (al^a^^) is equal to the square 
of the classical density and does not depend on k. 
Eq. (3.15) can be re-expressed as 



(«k«k) = 



2t[l + cos2 6*0(1 -7k)] cos So 



(3.16) 



from which one sees that the l/S'-corrected distribu- 
tion diverges like l/fik when approaching the condensate 
point at k = 0. More precisely, the momentum distribu- 
tion is singular and behaves like ^ sin9o/k for all values 
of the field —4 < fi/t < 4. The integral of this quantity 
over the whole Brillouin zone is convergent and yields the 
number of uncondensed particles, which is given by 



k#0 



+ ^Xl^'^°^^^"'^^k-^kWk) (3.17) 

k#0 

o(^k + WkVk) 



k^^O 
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-(S-> + 1/2 (LSW with corrected w.f) 
■ {S') + 1/2 (LSW without corrected w.f) 
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FIG. 1: (Color online) Hard-core boson density p as a function 
n/t. QMC results (symbols) are compared with classical and 
LSW calculations. The inset shows the relative deviation of 
the classical and spin- wave results from the numerically exact 
QMC estimates. 



Let us point out that one can also recover the conden- 
sate density using Eq. (B3) in Appendix B. Indeed, if k 
is replaced by zero in the right hand side of Eq. (B3), the 
expression can be rearranged to lead again to the con- 
densate density obtained previously using the definition 



IV. COMPARISON WITH QMC SIMULATIONS 

In this section we compare the large S approxima- 
tion of various quantities to exact QMC estimates ob- 
tained using the Stochastic Series Expansion (SSE) algo- 
rithm [24]. The simulations have been performed for the 
hard-core boson model Eq. (1.1) on square lattices Lx L 
with L = 8, 16, 24, 32 at temperatures low enough to 
get ground state estimates (/3/i oc L^), in particular in 
the dilute limit where the finite size gap scales as ~ 
with z = 2. Overall, we agree with the numerical results 
of Ref.[14] whenever we could compare. We have nev- 
ertheless included numerical results for the density, the 
condensate and the superfluid density to be able to dis- 
cuss their behaviour close to ^c, where the corrections 
to the harmonic ground state turn out to be crucial. In 
addition to these quantities that had already been dis- 
cussed in Ref.[14], this section also contains QMC results 
for the momentum distribution. 



A. Particle density 

Figure 1 is a plot of the hardcore boson density as a 
function of fi/t. The spin wave results for the density 
are plots of (S^) +1/2 calculated using the perturbed 
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0.25 



3. 



3 




0.15 



FIG. 2: (Color online) Logarithmic corrections to the den- 
sity as a function of distance from the critical chemical po- 
tential — ^c). QMC data (symbols) are described by a 
fit (dotted-dashed green line) of the form a |ln [e{^ — /ic)/4t]| 
with a ~ 0.034 and e ~ 0.04; the spin-wave calculation using 
the perturbed ground state (black line) captures the logarith- 
mic correction and yields a ~ 0.035 and e ~ 0.065. Classi- 
cal (blue dotted) and LSW results using the non-perturbed 
harmonic ground-state (magenta dashed) do not capture the 
logarithmic corrections. 



ground state (solid line) and using the harmonic ground 
state (dotted line). While both approaches yield signifi- 
cant corrections to the mean field density, the expectation 
value {S^) calculated in the harmonic ground state misses 
some terms of order 0(1), i.e. 1/5' corrections to the clas- 
sical results, as discussed in Sec. IIIB. This effect is best 
seen in the inset of Fig. (1), which shows the relative de- 
viation of the spin wave results from the QMC estimates. 
For small densities the relative deviation from the QMC 
result is as large as 40% if {S^-) is computed using the 
harmonic ground state. This deviation never exceeds 5% 
if the perturbed ground state is used. Furthermore, the 
computation in the non-perturbed harmonic ground state 
misses a very important feature of the dilute Bose gas 
limit. Indeed, logarithmic corrections have been shown 
to dominate the low-density limit [25-27] close to the crit- 
ical point /Ltc, and QMC data are indeed consistent with 
the behavior p ~ (/j, — /Zc) In (/j, — ^c) with = — 4i, 
as shown in Figure 2. The density computed using the 
perturbed ground state correctly captures the logarith- 
mic correction whereas {S^) computed in the harmonic 
ground-state does not. 

These results show without any ambiguity that the 
best way to estimate the density in the context of a semi- 
classical approximation is to deduce it from the derivative 
of the harmonic energy with respect to the chemical po- 
tential, or equivalently to deduce it from a calculation of 
the expectation value of in the ground state that in- 
cludes leading corrections beyond the harmonic approxi- 
mation. The density deduced from the expectation value 
of calculated in the harmonic ground state is much 
less accurate, and qualitatively wrong in the dilute limit. 




FIG. 3: (Color online) Superfluid density as a function of ^/t. 
Classical (dashed line), semi-classical (solid green line), and 
QMC (symbols) results are shown. Inset: superfluid fraction. 



^ ^ i i 




FIG. 4: (Color online) Condensate density as a function of 
fj,/t. classical (dashed line), semi-classical (solid green line), 
and QMC (symbols) results are shown. Inset: condensate 
fraction. 



B. Condensate density and superfluid density 

1. LSW and QMC results 

Figures 3 and 4 are plots of the condensate and su- 
perfluid densities as a function of The semi-classical 
results presented in the previous sections are in very good 
agreement with QMC results, as also discussed in Ref. 14. 
At the classical level, the condensate and the superfluid 
densities are equal. The effect of quantum fluctuations is 
to enhance the superfluidity and to deplete the conden- 
sate. QMC estimates for psf and po are obtained in the 
directed loop algorithm framework [24] using the winding 
number fluctuations [28] for p^i and the Green's function 
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estimate [29] for p^. Note that the latter suffers from 
larger statistical errors than p^i- 



2. Dilute Bose gas limit 

Building on the fact that the semiclassical results are 
very accurate, we analyze the extremely dilute limit for 
condensed and superfluid fractions using this approxi- 
mate framework with the help of semi-classical calcula- 
tions on finite square lattices of linear length L ~ 10^ 
down to very low particle density oi p = 10^^. Such a 
limit is simply impossible to access using QMC simula- 
tions, where the computational cost grows very fast, like 
~ L^, so that only systems with a linear size of the order 
of L ~ 10^ can be accessed. In both figures 3 and 4, these 
fractions are shown in the insets. While both fractions 
are the same at the classical level, the effect of quantum 
fluctuations is qualitatively different in the two cases. In 
the extreme dilute limit, they converge to 1 very differ- 
ently. Let us first consider the superfluid density. The 
superfluid fraction Psf/p is enhanced by quantum fluctu- 
ations with respect to the classical case, as seen in the 
inset of Fig. 3. More precisely, in the dilute limit, the 
semiclassical superfluid fraction (Fig. 5 left) tends to 1 
like: 



p.f/p=l-(cv)^ 



(4.1) 



with ( ~ 0.728, and an exponent v ~ 1.07 very close 
to one. Note that, at the classical level, both fractions 
(superfluid and condensate) tend to 1 like f = 1 — p. 

By contrast to the superfluid fraction, the condensed 
fraction is more affected by quantum fluctuations in the 
dilute limit. Indeed, it converges much more slowly to 
unity, as can be seen in the right panel of Fig. 5 and in 
the inset of Fig. 4. As first predicted by Schick in Ref. 25, 
logarithmic corrections of the form 



Po/p = 1 



|ln(CV)|- 



(4.2) 



are expected in the extreme dilute limit. Fig. 5 (right) 
shows semi-classical results for the very slow convergence 
of the condensate fraction to 1, with a fit to Eq. (4.2) with 
a ~ 0.86 and ^ ~ 0.68. Interestingly we observe that the 
effective distances ^ ^ ^ ^ 0.7. 



C. Momentum distribution 

We now turn to the momentum distribution 

N{k} = {ala^), (4.3) 

which can be efficiently computed using QMC simula- 
tions, following Ref. 29. Results for the half-filled case 
(/i = 0) are shown in Fig. 6 for k ^ along the line 
fcr — k,, = k in the first Brillouin zone. At the classical 




FIG. 5: (Color online) Semiclassical results (red symbols) for 
the superfluid (left) and condensate (right) fractions plotted 
versus the total density p. Blue lines are fits of the form 
Eq. (4.1) for the superfluid (left) and Eq. (4.2) for the con- 
densate (right). 



level (dashed line), the distribution does not depend on 
momentum and is equal to (pciass.-j2^ T\ve effect of spin 
wave fluctuations (solid line) is to introduce a momentum 
dependence which is singular near fc = and diverges like 
1/fc, as discussed in section IIIB. This redistribution of 
spectral weight is due to the fact that spin waves deplete 
the condensate at fc = 0. The semi-classical calculation 
of the momentum distribution reproduces the behavior of 
the QMC results for small fc. This is best seen in the right 
inset of the figure, which is a log-log plot of the distribu- 
tion showing the 1/fc dependence of the QMC results for 
small fc. Looking back at Eq. (3.16), the 1/fc divergence 
of iV(fc) is a consequence of the linear spectrum at small 
momentum fi^ ~ fc. Away from the condensation vector 
fc = 0, the agreement between semiclassical and QMC 
results is less good. The QMC estimate is consistent 
with a distribution that goes to zero when fc — >■ tt while, 
according to the semiclassical results, the distribution is 
only slightly renormalized downwards with respect to the 
classical constant value. 

The contribution at fc = is not shown on the main 
panel of Fig. 6 since it diverges with the system size like 
LP'. In the left inset however we show the QMC result 
for po — N{0)/L'^ plotted against 1/L. Using a quadratic 
fit, we extract the thermodynamic limit value of the con- 
densate density po — 0.188(2), in good agreement with 
the estimate 0.191(2) reported by Sandvik and Hamer in 
Ref. 30 (see also Table I). 



V. DENSITY SUM RULE 

Having obtained the first order corrections to the total 
density of particles, the density of condensed particles 
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FIG. 6: (Color online) Momentum distribution of hard-core 
bosons A'^(k) Eq (4.3) at half-filling (/i = 0) along the line 
= ky = k for fc 7^ 0. Classical (dashed line), semiclassical 
(full line), and QMC (different symbols for L = 16, 24, 32) 
results are shown together. Right inset: Log- log plot of the 
momentum distribution, which diverges as ~ 1/fc. Left inset: 
Finite size scaling of the condensate density from QMC as a 
function of 1/L. It is well accounted for by a quadratic fit 
(black line). 



and the density of uncondensed particles, it is natural to 
test the semiclassical approximation with respect to the 
following sum rule: 



1 V- t 



(5.1) 



The above equality follows directly from the conservation 
of the number of particles: the total number of hardcore 
bosons in the system is equal to the sum of the num- 
ber of condensed and uncondensed particles. Somewhat 
surprisingly, with the semiclassical expressions derived in 
the previous section, the sum rule of Eq. (5.1) is violated. 
Fig. 7 shows the average particle density obtained from 
Eq. (2.27) and a plot of the sum po + (Sk^o^il^k)/^ 
(sum of Eqs. (2.28) and (3.17)). In the inset, the vio- 
lation of the sum rule, defined as the relative difference 
between the two quantities, is shown in the entire filling 
range. It is smaller than 2% up to half-filling, and never 
exceeds 6.5% above half-filling. 

The origin of the violation of the sum rule is actually 
quite simple. In the spin-1/2 language, the sum rule relies 
on two identities: 



and 



1/2 = s+s; 



(5.2) 



(5.3) 




-4 



-2 



/i/t 



FIG. 7: (color online) Total density p (red diamonds) and p(,-\- 
X^icT^o '^k'^k (blue triangles) as a function of the chemical 
potential. The inset presents the violation to the sum rule. 



which lead to 



N ^ 



1 

iV2 



■'J 
kT^O 



E 



(5.4) 

since, according to the Matsubara-Matsuda transforma- 
tion, the local density is related to the z component of 
the spin by = S'f -I- 1/2. However, when the spin is 
larger than 1/2, the identity S'f + 1/2 = 5+5*^^ is no 
longer valid, and the expectation value of is no longer 
equal to that of 5*+ 5^ - 1/2. 

Remarkably enough, in spite of that, the sum rule is 
almost satisfied, especially below half-filling. This pro- 
vides some additional confidence in the accuracy of the 
semi-classical approach. As a further test, we discuss in 
the next section the effect of higher order corrections on 
the ground state energy. 



VI. SECOND ORDER SPIN WAVE THEORY 
A. Energy 

The convergence of the expansion of spin operators in 
terms of Holstein Primakoff bosons Eq. (2.8) can be a 
cause of concern in the case S = 1/2. In this section, we 
compute the 1/5'^ correction to the expectation value of 
several observables. 

The ground state energy of the Hamiltonian to 'HS^^ 
gives the energy of the original problem to order XjS. 
The 1/5*2 correction comes from both H^^) ^(4) j^. 
is obtained by treating 'H'^-' up to second order in per- 
turbation theory, and "H^^^ to first order, as pointed out 
in another context by Zhitomirsky and Nikuni[31]. The 
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contribution of "H*^^^ is simply given by E^^^ = ('H^^^) /N, 
where the expectation value is calculated in the harmonic 
ground state. Using Wick's theorem, this contribution is 
given by: 



= lim — - [cos^ e(2m ~6)(n~ A) 

+ (2m + 5)(A + n) - sin^ e{Am^ + 



(6.1) 



where m,n,S and A are defined by: 



^^"k^'k = A=^E 



(6.2) 

where i and j are nearest neighbors. 

Being odd in the number of bosonic operators, 'H'"^-' 
contributes only at second order in non degenerate per- 
turbation theory: 



^(3)^1 ^ l(e|H(3)|0)P 



N 



\e)r^\0) 



|0> 



E, 



o 



(6.3) 



where E''^-' is the energy per site. The effect of 'H^^-' on 
the Bogoliubov vacuum is to create either single magnon 
excited states or three magnon excited states. We treat 
these two cases independently and write the H*^'^^ contri- 
bution to the energy as E^^'' — e[^ 
magnon component of H^'^^O) is: 



4^. The single 



2< sin cos 6* ■r-^ n rr, 2 / 2 n \ 

V V 

(6.4) 



which leads to: 



E 



(3) 



lim —r- sm cos h 



E [2wk + 7k(t'k - ^iki'k)] 



1 

- -4t7 cos^ 6*0 hm (2m + n - A)^ 
25^2 r^o 



(3) 

In the thermodynamic limit E^^ is dominated by the 
excited states in which the three magnons all have differ- 
ent momenta. The three-magnon component of 'H*''^'|0) 
that fulfills this condition is given by: 



.5) 



2t sin 9 cos 9 

~sVM7W 



E'f(q,k)|lklql-k-q) 



k.q (6.6) 

f(q: k) = 7q("k-|-qWkWq - "k^q^k+q) 

where the sum ^ is such that the three momenta 
k, q and — k — q are all different. The three-magnon 



1 order SW 
2"" order SW 
• • QMC 




FIG. 8: (Color online) Plot of the 1°* order spin wave energy, 
of the 2°"* order spin wave energy and of the QMC energy 
measured with respect to the mean field energy. 



contribution to E^^^f takes the form 
~2t^ 1 



E. 



(3) ^ 
3m. 



53 



sin^ 9q cos^ 9o 



1 



■ lim — — — 
r^o 3! ^ 17k 

k,q 



F'(k,q) 



■7) 



fin + illi 



'q T- ^'k+q 

where F(k,q) is defined by: 

F(k,q) = f(k,q)+f(-k-q,q)+f(q,k)+f(q,-k-q) 
+f(-k-q,k) + f(k,-k-q). 

(6.8) 

The spin- wave approximation of the energy per site to 
order 0{1/S^) 



E = £ + £;(2) + e[ 



(3) 



E. 



(3) 



(6.9) 



is plotted in Fig. (8), together with the 1"* order SWT 
and QMC results, as a function of fj,/t. All energies are 
measured with respect to the classical energy £. The 
difference with the classical energy is monotonously in- 
creasing (in absolute value) from the dilute limit up to 
half-filling, where the correction due to quantum fluctu- 
ations is the most important. Clearly, the 1/S correction 
captures most the quantum correction, but the inclusion 
of 1/52 corrections leads to a significantly better agree- 
ment with QMC results. This systematic improvement 
upon including higher order 1/S corrections gives addi- 
tional support to the semi-classical expansion. 



B. Other observables 

Superfluid and condensate densities, as well as the 
compressibility k = dp/dfi can also be calculated at or- 
der 1/5*2. The calculation is straightforward but cumber- 
some, and for simplicity it is not reproduced here. The 
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estimates that we have obtained at half-filhng are listed 
in Table I, together with QMC results from the present 
work as well as from Ref. 30. For the ground-state energy 
and the compressibility, the agreement increases system- 
atically from 1/5*° (classical) to l/S"^. For the energy, 
the relative error is ^ 9% for the classical estimate, less 
than 1.3% including the 1/S correction, and ~ 0.2% in- 
cluding the l/S"^ correction. For the condensate density, 
the result including the 1/S' correction is already within 
the error bars of QMC, and it is not clear whether in- 
cluding \/S'^ correction leads to any improvement. For 
the superfluid density, both the results including correc- 
tions up to order 1/S and lie outside the error bars 
of QMC, and the result up to order 1/S appears to be 
better than the result up to order 1/S'^. In any case, 
the improvement over the classical result is clear for all 
quantities. 







Psf 


PQ 




Classical 


-1 


0.25 


0.25 


0.125 


1/S sw 


-1.08382 


0.27095 


0.18904 


0.1075 


i/s^ sw 


-1.09539 


0.27198 


0.19127 


0.1053 


QMC (Ref. 30) 


-1.097648(4) 


0.2696(2) 


0.191(2) 


0.1048(1) 


QMC (this work) 


-1.09764(1) 


0.2697(2) 


0.188(2) 


0.1048(1) 



TABLE I: Ground-state estimates at half-filling (/i = 0) for 
the energy per site eo, the superfluid density psf, the conden- 
sate density po, and the compressibility The three first 
lines are analytical results from classical and spin-waves at 
first (1/5 SW) and second order (l/S^ SW). Below are shown 
QMC estimates from SSE simulations obtained by Sandvik in 
Ref. 30 and in this work. 



VII. CONCLUSION 

The semi-classical approach to hard-core bosons on a 
lattice, which is based on a large S approximation to 
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the Matsubara-Matsuda spin-1/2 version of the Hamil- 
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What is the correct way to get the exact 1/S correction 
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Appendix A: Computation of the superfluid density 

1. Classical value 

The superfluid density can be obtained by imposing a 
phase gradient $ri+e ~ — to the system (e being 
the unit vector along the x or y axis of the lattice). At 
the classical level, this leads to the following energy cost 
per site 



E{^) - E{0) = t (sin6lo) if^ + 0(/). 



(Al) 



Using the analogy introduced by Fisher, Barber and Jas- 
trow in Ref. 32 where the kinetic energy density of a 
superflow (density psf and velocity Vgi) in one direction 
is SE{vs() = lm*ps{ 
gives for the superfluid density 



- (wgf)^, with Vsf — {h/m*)ip. This 



2h^ 



^sf, 



(A2) 



where Tgf = d'^E{ip)/dip^\^^^ is the helicity modulus, 

the effective mass is given by 2m*//i^ = l/{2t), and the 
factor of 2 comes from the fact that the twist (p has been 
introduced in both directions. Finally we get at the clas- 
sical level 



sin^ 9n 



Psf 



(A3) 



Interestingly, we remark that condensate and superfluid 
densities are equal at this level of approximation 



SW corrections 



In order to evaluate the SW corrections to the super- 
fluid fraction, the phase gradient can be introduced di- 
rectly on the bosonic operators 



at. e^*-» 



Or. ^ a,^ e-'*'. , (A4) 
which in term of equivalent spin operators translates into 
S^^ S^^ cos - SI sin $r. 

SI S^^ sin + SI COS $r. ■ (A5) 

Therefore, the rotation (2.7) becomes 

S^^ = (cos eS^ + sin eS^^ COS - S^ sin $r. 
5*^^^ = (cos 9 SI + sin OS^^ sin $r. + S^ cos $r. 

(A6) 



-sines^.+coses:! 



In the new rotated frame, at the linear SW approxima- 
tion the XY Hamiltonian now reads 

H^^\ip) = 2^[Ak((/5) (a|,ak + a^_ka-k 



(A7) 



14 



with 



M^) = "2 cos 95 [(1 + cos^ 0)7k - 4] , (A8) 



The terms involved in {S^ Sj ) up to order 0{S) are: 



and 



^kCv) = 2 cos(p(sin^6')7k, 



(A9) 



where we used S—1/2 and where 9 is fixed by the equa- 
tion COS0 = fi / {4t cos (j)) imposed by the minimization 
of the classical energy. At order ip"^ we have of course 
cost/? ~ 1 — (p'^/2. Similarly the condition on 9 yields 
cos^ 9 ~ cos^ 9o (1 + (/3^), and sin^ 9 ~ sin^ 6'o — cos^ 9oip^. 
Therefore we have up to the order ip'^: 



^k(<p) = ^k(O) - ^ (Ak(0) + t7k cos^ 9o) (AlO) 



and 



Sk(</5) = Sk(0) - ^ (Bk(0) + t7k cos^ 9o) . (All) 



(SrS-) = ^cos'9{0\{h, + bl){b,+b])\0) 



S 

sin^ 61 - S-sin^ 9{0\blb, + b]bj\0) 
-'^mh-bl)ib,-b]m 



- cos a sm 



0— ((0|(6. + 61)10^) +h.c. 



((O|(6,+fo])|05 



5) +h.c. 



-5 cos ( 



(B2) 

Injecting this result in Eq. (3.14) and making use of the 
definitions of the inverse Fourier transforms of the bi op- 
erators, Eq. (2.13), we obtain: 



Writing ^k(O) Ak and Bk(0) Sk, the 1/S correction 
to the GS energy in the presence of a small twist reads 



t7kCos 9o 



V Ak + SkJ' 



[a^ai,) - hm {| cos^ 0(O|(6_k + 6l)(6k + folk)|0) 
+8"^ sin^ 6l7V(5k,o 

-5sin2 0,5k,o(O| J2 bhie^'^''' + e-^^^OlO) 



Thus the superfluid density is given at 1/S order by the 
following expression 



Psf 



sin^ 9q 
4 

1 

mt 



i7k cos^ 9q 



V^k + BkJ' 



(A12) 



This expression and that of Ref. 14 are strictly equiv- 
alent only at n = —4:t (low density limit) and at = 
(half fining). In the range -4 < n/t < and < fi/t < 4 
they differ by a term which is of order [33] . 



Appendix B: momentum distribution 



We start by expressing S^Sj in the rotated frame: 



5*+ 5*7 cos2 9Sfsf + sin^ 9Sf Sf + sf S^' 



■ cos f sm ( 



(Bl) 



-cosflS'f Si J — sin 9 Sf 6ij. 



--(0|(6_k-6|,)(6k-&lk)|0> 



-I- cos 9 sin ( 



NS 



+ cos sm ( 



V2 
TVS' 



5k,o ((O|(6_k + 6[)|05)+h.c. 



0sin0^^Jk,o ((0|(6k + &lk)l'^^)+h.c.) 



-t-S'cos6'|. 



(B3) 

Hence, the number of particles at momentum k 7^ is 
given by 



(OkOk) 



hm {| cos^ 0(O|(5_k + &D(&k + 6lk)|0> 
mb-i,-bl){by,-bl^J\0)+Scos9} 



-(l + COS^o)' 

[(1 + cos^ 9o)vl + MkWk sin^ 6*0] /2. 



(B4) 



